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ABSTRACT 

This study concerns the fast and accurate solution of multilevel non-LTE radiation transfer problems. We propose and evaluate an 
alternative iterative scheme to the classical MALI method. Our study is indeed based on the application of Broyden's method for the 
solution of nonlinear systems of equations. Comparative tests, in ID plane-parallel geometry, between the popular MALI method and 
our alternative method are discussed. The Broyden method is typically 4.5 times faster than MALI. It makes it also fairly competitive 
with Gauss-Seidel and Successive Over-Relaxation methods developed after MALI. 
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1. Introduction 

The solution of the non-LTE multilevel-atom radiative transfer 
problem is a classical one in astrophysics. Indeed, the assump- 
tion of non-LTE implies, consistently with the departure of the 
source functions from Planck functions, that the population den- 
sity of the atomic or molecular levels considered depart from 
what can be derived at LTE, in a straightforward manner, using 
Saha and Boltzmann relations (see e.g., Mihalas 1978). 

In the non-LTE case, one has on the contrary to solve simul- 
taneously and self-consistently for a set of A^x equations of ra- 
diative transfer together with A'l equations of statistical equilib- 
rium (hereafter ESE) describing the detailed balanced between 
excitation and de-excitation processes between every atomic or 
molecular levels. Since absorption and stimulated emission ra- 
diative rates depends explicitely on the radiation field, which 
itself depends on the level populations, this problem is intrin- 
sically a search for the solution of coupled nonlinear equations. 

Since the beginning of numerical radiative transfer in the 
late 60's, the two most popular methods used for tackling this 
problem have been the complete linearization method of Auer & 
Mihalas (1969) and the Accelerated A-Iteration based scheme 
called MALI (Rybicki & Hummer 1991). Despite their apparent 
differences, they have however in common the fact that basically, 
one is conducted to deal with linearized equations. An interest- 
ing comparative study of these two approaches have been made 
by Socas-Navarro & TrujUlo Bueno (1997). 

In this study, we investigate on the use of a quasi-Newton 
numerical method for the solution of the nonlinear ESE. Our 
choice went to Broyden's method (1965) whose elements will be 
presented in §2. To the best of our knowledge, Koesterke et al. 
(1992) were the first to bring this numerical scheme into the field 
of radiation transfer. Their study was presented in the context of 
the modelling of spherically expanding atmospheres of hot and 
massive Wolf-Rayet stars. Broyden's method was more recently 



invoked in the context of the coupled-escape probability method 
(Elitzur & Asensio Ramos 2006). 

Besides from the required algebra and mention to caveats 
related to the implementation of the method, it remains how- 
ever difficult to figure out from Koesterke et al. (1992) the ac- 
tual performances of such an approach. A comparison with an- 
other method have also been barely evoked by the authors, who 
mentioned however a significant speed-up provided by Broyden 
algorithm for large A^l atomic models. In particular, being con- 
temporary with the publication of Rybicki & Hummer (1991), 
it is a pity that no comparison with the MALI method could be 
made yet. Such an evaluation is the scope of the present work. 

2. The numerical scheme 

For a A^£-level atomic model, the ESE will in general write as a 
set of elementary equations: 

^[n;Aij - {njBji - niBij)Jij] 

j<' 

j>i 

+ Yj^niCij-njCji)^0 (1) 

i*i 

where the A,^ and By stand respectively for the spontaneous 
emission, and the absorption and stimulated emission rates, n, 
represents the population density for each energy level, and Jij 
is the scattering integral for each radiatively allowed transition 
we shall consider. Besides the radiative processes, the C,j are 
coUisional excitation and de-excitation rates. In general, these 
rates depend on the electronic density so that, if the latter is not 
known a priori, terms like n,C,y are nonlinear in the population 
densities. Hereafter we shall consider only cases for which the 
coUisional rates are known a priori. 
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Fig. 1. Typical relative error, (thin), convergence error, Ck 
(thick) and (dotted) respectively for MALI (dash) and 
Broyden (full) vs. the number of iteration, for both scheme. 

The scattering integral entering the ESE is formally written 



Jij = Aij[Sij]^ 



(2) 



where, assuming complete redistribution in frequency, the 
source function is defined as: 



HiBij 



(3) 



A large system of Eqs. ([T]) is homogeneous so, practically, 
one of the equations have to be replaced by a constraint equation 
like, for instance, a conservation equation of the form: 



Nl 

>=1 



(4) 



2.1. Broyden's algorithm 

The system of Eqs. ([T]) and Q can be reformulated by defining a 
function F acting on the set of n*'"* = (Wj'^', \ njy^) where 
T is a discrete depth along the opacity scale used to sample our 
slab or atmosphere. F is defined such as: 



Ff'- Z [n'pA.j-UpB.-n^pB.yj^j] 
- 2 [nf - ir^pBij - n^Bj^Jij] 



+ 2 («r C,7 - n'l'Cji) 
for / ^ Nh and, if / = A^l 

Nl 



(5) 



(6) 



Computations of F requires repeated evaluations of scattering 
integrals Jjj - Eq. (2) - that we perform here using the well- 
know short-characteristics method with monotonic parabolic in- 
terpolation introduced by Auer & Paletou (1994). 



O 



-10 - 
1.0 



1 












\ 


* 




\ 


* 






\ t 









2.0 

log(time [s]) 



Fig. 2. The convergence error, Ct, for MALI (dashed) and 
Broyden (full) are displayed vs. computing time for a 5-level H 
atom and, respectively, 5 (thin) and 8 (thick) points per decade 
in optical depth. 



In that frame, and using the Sherman-Morrison formula (see 
e.g.. Press et al. 1992) which provides an analytical formula for 
the direct computation of the inverse of the Broyden matrix, our 
algorithm consists in the following steps. We choose an initial 
vector no, at every t depth, and an initial Broyden matrix Bq; 
then we compute B^K The iterative scheme is such that: 



6nk = -B^ F(nk) , 

then update 

n(k+i) - nk + 6nk , 

then compute 

6Fk = Fin^k+i)) - F(nk) , 

and finally, update 



-I , (Snk-B,'SFk)6nlB-^ 



B(k+\) - Bk + 



{SnlB-^')6Fk 



(7) 



(8) 



(9) 



(10) 



while IIFII2 > B. Practically, we have chosen e = 10 ^, which 
guarantees that we have reached a fully converged state (see also 
Fig-Q. 

2.2. Initialization 

The proper initialization of the Broyden scheme is a critical is- 
sue. We employed the following method which was tested as 
suitable, both from the standpoint of an adequate start and from 
the one of an acceptable computing time. 

Before starting the iterative scheme, we assume LTE popu- 
lations for our model-atom. In such a way, the grand function F 
defined by Eqs. (jSj) and (j6]l can be fully evaluated. Then, we com- 
pute an initial Jacobian, Bq, using the finite difference scheme 
fdjac described by Press et al. (1992). 

In the following figures relative to the timing properties of 
Broyden's method vs. MALI, we shall always include the spe- 
cific time necessary for the evaluation of Bq. 
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Fig. 3. Respective computing times vs. number of levels of the 
H-atom model, for for MALI (dashed) and Broyden (full). The 
dash-dotted curve corresponds to the time required for the eval- 
uation of the initial Jacobian. 



3. Comparison of Broyden vs. MALI 

We adopted the popular "flavor" of the MALI method, us- 
ing a diagonal approximate operator, as described by Rybicki 
& Hummer (1991), and without acceleration of convergence 
schemes. 

In order to compare the properties of the Broyden scheme vs. 
MALI, we adopted a ID semi-infinite plane parallel slab model 
of Tmax - lO'", discrctizcd in various number of points per 
decade in optical depth, using also 3 to 10 energy levels H-atom, 
inspired by the classic benchmark proposed by Avrett (1968; see 
also Leger & Paletou 2007). As in the latter, the slab tempera- 
ture was fixed at 5000 K and the collisional rates set at 10^ s"'. 
We also adopted the definitions initially proposed by Auer et al. 
(1994) for the relative error, from an iteration {k) to another: 



n(k) 



and for the "convergence error": 



'(oo) 



(11) 



(12) 



where n^ca) is the "fully converged" solution obtained, for a given 
method and model, after a large number of iterations. We also 
introduce the quantity 



Fk = Fwlb 



(13) 



i.e., the Euclidian norm of F, the function defined by Eqs. (|5]l and 



\6\. Note that f',^'' for the MALI method is defined by a modi- 



ed Eq. (j5]l following the preconditioning strategy proposed by 
Rybicky & Hummer (1991). 

3.1. Convergence 

In Fig.[T] we display the rates of convergence of the Broyden and 
MALI methods, respectively. The convergence error C/t for each 
method have been computed using population densities obtained 
once Fk < 10"^ in both cases. The case used here is a 5-level H 



atom with a ID slab discretized by a 5 points per decade grid 
in ojjtical depth. It is worth noting that, in order to reach Fit < 
10^' and a well-converged solution though, one should iterate 
up to reaching Rii values as small as 10"'" typically. In terms 
of number of iterations, Broyden typically beats MALI of more 
than one order of magnitude. However, the quite distinct nature 
of each method makes such a comparison incomplete. Hereafter, 
we carry on such an analysis but we shall compare respective 
computing times. 

3.2. Sensitivity to tiie spatial (optical depth) refinement 

In Fig.|2] we turn to an analysis of the respective timing prop- 
erties of Broyden and MALI. It is shown that Broyden, again, 
always beats MALI by a typical factor of the order of about 4- 
5 in time. It is also the case when the spatial grid refinement is 
increased from 5 to 8 points per decade, for instance. It is impor- 
tant to note that timings given for Broyden include the computa- 
tion of the initial matrix Z?o. This is why the rates of convergence 
displayed for the Broyden method do not start at f = in Fig.|2] 

3.3. Sensitivity to the number of transitions 

A next important point to investigate relates to the advantage of 
Broyden against MALI when an increasing number of atomic 
transitions is considered. Again, as demonstrated in Fig.|3] the 
full Broyden iterative process is always significantly faster than 
MALI. In general, the gain due to the Broyden method is of the 
order of 4-5 in total computing time. This is less than the the gain 
of the order of 8 already reported by Koesterke et al. (1992), al- 
though their method of reference was presumably different from 
MALI. 



3.4. Discussion 

We are aware that the MALI method can be speeded-up by accel- 
eration of convergence schemes (see e.g., Auer 1991). However 
the most significant improvements in the field of iterative meth- 
ods for the non-LTE radiative transfer problem were brought by 
the introduction of the Gauss-Seidel (GS) and Successive Over- 
Relaxation (SOR) methods (Trujillo Bueno & Fabiani Bendicho 
1995). It was already shown, for instance, that SOR always beats 
both Jacobi (i.e., accelerated A-iteration with the diagonal of the 
full A operator as an approximate operator) and GS, even when 
Ng acceleration of convergence is applied. 

Beyond the fact that Broyden is significantly faster than 
MALI, we can also add that Broyden is as competitive as the 
SOR method, according to Paletou & Leger (2007; see their 
Table 1 where comparable timing and the corresponding itera- 
tion numbers were given for MALI, as we used it in the present 
study, GS and SOR). 

The Broyden method is also potentially more advantageous 
than MALI and GS/SOR, because of its intrinsic capability to 
deal with the self-consistent evaluation of the electron density 
in a multilevel non-LTE problem, if necessary - a problem for 
which MALI needs to plug-in a Newton-Raphson scheme to it, 
as proposed by Heinzel (1995) and Paletou (1995). 

Another important point is that, as indicated in our Fig. 3, a 
great deal of time of our Broyden code is spent in the computa- 
tion of the initial Jacobian, a task which can be performed with 
great advantage using parallel computing. The inner structure of 
the fdjac routine permits, indeed, parallelization with a high 
scalability. 
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As a final comment, it is also important to consider that 
Broyden's method can be easily implemented in already exist- 
ing codes, without the need of modifying the formal solver, un- 
like with GS/SOR methods. 

4. Conclusion 

We propose an alternative method for the solution of the 
non-LTE multilevel radiation transfer problem. It is based on 
Broyden's method for the solution of nonlinear systems of equa- 
tions. The method is easy to implement and it is about of factor 
of 4.5 times faster than the well-known MALI method. Another 
advantage is that it does not require any modification of usual 
formal solvers, as it is the case for GS-SOR methods developed 
after MALI. It is also potentially very well- suited for parallel 
computing. Further tests will include the self-consistent treat- 
ment of the ionization balance, usually treated together with 
MALI with a Newton-Raphson scheme. In a next step, we shall 
consider more demanding models such as H2O, for instance. 
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